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Abstract 

To study the possibility of a fluid-fluid phase transition, we analyze a three- 
dimensional soft-core isotropic potential for a one-component system. We 
utilize two independent numerical approaches, (i) integral equation in the 
hypernetted-chain approximation and (ii) molecular dynamics simulations, 
and find for both approaches a fluid-fluid phase transition as well as the 
conventional gas-liquid critical point. We also study the possible existence 
of a triple point in the supercooled fluid phase at which three phases — gas, 

high-density fluid, and low-density fluid — coexist. 
PACS numbers: 61.20.Gy, 61.25.Em, 65.70.+y, 64.70.Ja 
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Recent experiments on phosphorus j|] show that a fluid-fluid phase transition is present 
at a temperature above the gas-liquid critical temperature. The transition is between two 
stable fluids at different densities. This sort of transition between low-density (LD) and high- 
density (HD) fluids has been proposed for many different materials, including supercooled 
H 2 0, AI2O3-Y2O3, Si0 2 , Ge0 2 , C, S, Ga, Se, Te, I 2 , Cs, Bi and Si 0. Experimental and 
theoretical studies support this hypothesis, and some suggest (e.g. for H 2 and C) the 
presence of a fluid-fluid critical point (CP) [0,|H3l. 

Theoretical and numerical studies take into account, as accurately as possible, the 
anisotropies of the molecular interactions, depending on the relative molecular orientations 
0. An alternative approach to the problem consists of investigating simple isotropic poten- 
tials which might be able to give rise to the LD-HD fluid transition. This makes it possible 
to find out which are the ingredients of the inter-particle interaction related to LD-HD fluid 
transition. 

Here we study an isotropic potential in three dimensions (3D) using two methods: (i) 
integral equations in the hypernetted-chain (HNC) approximation and (ii) molecular dy- 
namics (MD) simulations. These approaches agree qualitatively, showing a phase diagram 
with a gas-liquid CP and a LD-HD fluid transition with a second CP. The results suggest 
the presence of a gas-LD-HD triple point. The MD simulations reveal that both transitions 
are in the supercooled fluid phase, i.e. where the fluid phase is metastable with respect to 
the solid. 

Stell and Hemmer || proposed the possibility that an isotropic potential with a re- 
gion of negative curvature in the repulsive core ( "core-softened potential" ) might have two 
transitions. Recently it was shown, using analytic solutions in ID and MD in 2D, that a 
Stell-Hemmer type of interaction gives rise to water-like fluid anomalies that can be related 
to the existence of two different local structures in the fluid phase ||. Furthermore, the 
behavior of spherical particles, interacting in 3D through an isotropic potential with a hard- 
core and a linear repulsive shoulder, has been investigated through Monte Carlo. [|TtJ. It 
was also shown that the addition of an infinite range van der Waals type attractive term 
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may give rise to a second CP. The attraction, however, was not explicitly included in the 
simulation and was taken into account in a mean field type scheme (jTD) . 

One feature of LD-HD transition is that there is a competition between an expanded 
structure and a compact structure. In water, for example, the former is due to hydrogen- 
bond formation and is preferred at low pressures and low temperatures, while the latter 
is favored at high pressures and high temperatures. This suggests a potential with two 
characteristic radii ||. Here we consider a 3D potential (Fig. [I], inset) with a hard-core 
radius a and a soft-core radius b. Depending on the value of pressure and temperature, the 
effective hard-core will be a or b. The soft-core is given by a repulsive square shoulder with 
height Ae. The attractive part of the potential is given by a square well, with depth e and 
an interaction cut-off at distance c. Choosing a as length unit and e as energy unit, the 
potential has three free parameters: A (we set e = 1) and ratios b/a, c/a. Their possible 
combinations are too many to be investigated directly by means of MD. Therefore, it is 
advisable to use a faster, way of estimating the phase diagram. In this way it is possible to 
have a first indication of the phase diagram dependence on the potential's parameters. This 
allows one to perform MD simulation in an systematic fashion. 

Integral equations ]TJJ for the radial distribution function may provide a convenient 



tool for estimating the phase stability of the fluid. In particular, the HNC equation fll 
has often been used for relating phase diagrams and potentials |12|]. For this theory there 
exists a region in the density-temperature (p-T) plane in which no solution can be found. 
When this region is approached from high temperatures, there is a strong increase of the 
isothermal compressibility Kt- Since there is no genuine divergence in K?, this region 
cannot be identified in a rigorous way with the spinodal decomposition region of the fluid. 
However, for a large variety of potentials, its shape resembles qualitatively the fluid spinodal 



| T2|| . This allows us to have a first hint of the topology of the fluid instability region. 

The instability line, below which no solution of the HNC can be found, is shown in Figs. [I], 
|2|. For b/a = 1.4, c/a = 1.7 and A = —0.5, the potential considered in Ref. ||, the shape of 
the instability line is similar to that of a standard gas-liquid spinodals and does not change 



significantly by varying the value of A. However, for other combinations of the parameters, 
the shape can change dramatically (Fig. |2[). For example, for b/a = 2 and c/a = 2.2, we 
have a bifurcated instability line at intermediate values of A. The HNC instability line is 
not a true spinodal, but its bifurcation, with two distinct maxima, supports the existence 
of two CPs JT3[]. Note that the bifurcation in the instability line disappears by increasing 
A (as shown in Fig. |2|.A) or by changing only one of the other parameters. Therefore, the 
phase diagram depends drastically on the delicate interplay between the parameters. 

Based on the HNC results, we performed extensive MD simulations for a system of 
particles interacting through the potential with b/a = 2.0, c/a = 2.2 and A = 0.5 At 
high density the stable phase is solid, but the MD allows us to study the supercooled fluid 
phase before the crystallization occurs. The life-time of this metastable phase depends on 
(p-T). To be sure that our estimates are done in the fluid phase, we study the structure 
factor S(Q, t) = [p*At)pq(t)]/ N where p^ = J2j exp(iQ ■ rj(t)), N is the number of particles, 
Q is the wave vector [TJj and fj(t) is the position of particle j at time t. The growth of 
S(Q,t) for wave vectors with Q — > indicates a phase separation. Furthermore, for finite 
Q, it marks the onset of the crystal phase. Our simulations show the occurrence of a phase 
separation without a crystal formation (time interval labeled with B in Fig. |3|), followed by 
the onset of crystal phase (C in Fig. [3]) marked by a dramatic increase of S(Q,t) for Q — 6 
and 12 (in a -1 units) ||19|| . Our calculations are done in the (metastable) fluid phase before 
the crystallization. 

The study of the static structure factor S(Q) = [(S(Q,t))], where (•) is the average 
over Q vectors, and [•] is the average over a time interval, done on different time intervals 
(A, B and C in Fig. |3|), confirms the previous analysis. To emphasize the phase separation 
occurring in the fluid phase, we show in Fig. || the snapshot corresponding to the larger peak 
in time-interval B in Fig. [|a. The histograms of the number of particles in the system for 
the three planar projections of the snapshot show a clear separation in density. Therefore, 
the phase separation marked by the increasing of S(Q) at low values of Q is associated with 
the coexistence of LD and HD fluids. 
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The MD phase diagram is shown in Fig. [5[ At low densities and high temperature, the 
system is in the gas phase. Increasing the density or decreasing the temperature, the system 
undergoes a first order gas-solid transition. In the supercooled fluid phase, we observe two 
regions with negatively-sloped isotherms. The first region is for temperatures ksT/e < 0.61 
in the low density regime. The second region is for temperatures ksT/t < 0.645 in the 
high density regime. These regions correspond to thermodynamically unstable states and 
their borders are the spinodal lines. In the region of low density the coexisting fluids are 
the gas and the liquid. In the region of high density the analysis presented above allows 
one to conclude that the coexisting phases are a LD and a HD fluid. The points at highest 
pressure on the spinodal lines correspond to the CPs. In both HNC and MD approaches 
the merging of instability regions below the critical temperatures suggests the merging of 
coexisting lines, i.e. the existence of a gas-liquid-liquid triple point. 

Within our precision we find that the gas-liquid critical point is at ksT c i/e = 0.605 ± 
0.005, a 3 P c i/e = 0.018 ± 0.002 and a 3 p cl /m = 0.115 ± 0.015 corresponding to a packing 
fraction r] cl = 0.06 ± 0.01 (rj = ira 3 p/6). Note that T cl and p c \ are in good agreement 
with the HNC results. The second instability region occurs at packing fraction rj > 0.08, 
deep in the solid phase, where the metastable fluids phases have short lifetimes. Within 
our precision the critical temperature is ksT^/e = 0.650 ± 0.005, the critical pressure is 
a 3 P C 2/e = 0.08 ± 0.01 and the critical packing fraction is r] c2 = 0.16 ± 0.01. Note that T c2 is 
consistent with the HNC result, while p c2 is smaller. 

Since the MD calculations show that T c2 > T cl , the above set of parameters is not able 
to reproduce a fluid-fluid critical point. Nevertheless the HNC analysis suggests that other 
sets of parameters can lead to T c2 < T c \ as shown in Fig. 0.b. Preliminary calculation on 
a second set of parameters with b/a = 2.2, c/a = 2.4 and A = 0.5 show that both CPs are 
lowered in temperature and pressure. 

In conclusion, we studied an isotropic soft-core potential with a repulsive shoulder and 
an attractive well by means of integral equations in the hypernetted-chain approximation 
(HNC) and molecular dynamics (MD) simulations. Both approaches show, in the super- 
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cooled fluid phase, a gas-liquid critical point (at low density, low pressure, low temperature) 
and a transition (at high density, high pressure, high temperature) between a low-density 
(LD) and a high-density (HD) fluid with a second critical point. This second transition 
resembles the fluid-fluid transition discovered recently in experiments with phosphorus [|TJ. 
However, in phosphorus the fluid-fluid transition occurs between stable phases, while in our 
case it is in the supercooled phase as in H2O. Furthermore, both numerical approaches used 
here suggest the existence of a gas-liquid-liquid triple point. Therefore, the potential pro- 
posed here, despite its simplicity, is able to reproduce many intriguing features of materials 
with a second fluid-fluid transition, leading to a new way to explore the existence of a second 
critical point. 

We wish to thank P.V. Giaquinta, G. Pellicane, A. Scala, F. Sciortino, F. W. Starr, 
for helpful suggestions and for interesting and stimulating discussions. We thank NSF for 
support and CNR grant N.203. 15.8/204.4607 (Italy) for partial support. 
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FIGURES 

FIG. 1. Instability line in (p-T) plane of the HNC equation for b/a = 1.4, c/a = 1.7 e = 1 and 
A given by the label of the curves. Inset: general form of the soft-core potential. 

FIG. 2. As in Fig.l with (A) b/a = 2.0, c/a = 2.2 and (B) b/a = 2.5, c/a = 3.0 and A given by 
the label of the curves. 

FIG. 3. Structure factor versus time for (a) wave vectors with Q ~ 1 and (b) Q ~ 12 (in a -1 
units) at density a^p/m = 0.27 and ksT/e = 0.62 up to 3 x 10 4 UPP; (c) The static structure 
factor S{Q) averaged over different time intervals: labels (A,B,C) stands for the time intervals 
shown in panels (a,b). For A (B) we averaged over a low (high) density fluctuation time interval; 
for C over an interval after the crystal seed formation. The curves for S(Q) are shifted for the sake 
of comparison. 

FIG. 4. Histograms of number of particles as function of one coordinate for the snapshot 
corresponding to the large value of S(Q = l,t) in time interval B of previous figure. The corre- 
sponding projections of the snapshot are also shown. Histogram and snapshots are shifted for the 
sake of clarity. The horizontal line shows the average density. 

FIG. 5. Pressure-density isotherms from the MD simulation for the potential with b/a = 2.0, 
c/a = 2.2 and A = 0.5. Full points represent stable fluid. Open points represent supercooled fluid 
[15]. The inset emphasizes the instability region with negatively-sloped isotherms (and a van der 



Waals loops at ksT/e = 0.6) associated with the gas-liquid CP. It is possible to see a high density 
region with the van der Waals loops associated with the second CP. The lines connecting the points 
are just guides for the eyes. Temperature is measured in e/ks- 
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